 // ==+==+==+==+==+==+==+==+==+==+==+==+==//
 // Modelo de Simulacion Sardina  (MOSAR) 
 // VERSION 1.0     2011-10-25
 //  LA Cubillos - EPOMAR
 // Universidad de Concepcion
 // Demandante: Subsecretaria de Pesca
 //======================================
 //Bitacora LOG
 //

DATA_SECTION
	!! cout<<"MOSAR detecta la plataforma "<<PLATFORM<< endl;

	init_adstring DataFile;
	init_adstring ControlFile;

	!! BaseFileName=stripExtension(ControlFile);
	!! cout<<BaseFileName<<endl;
	!! ReportFileName = BaseFileName + adstring(".rep");

	!! ad_comm::change_datafile_name(DataFile);

	int SimFlag;
	int rseed;
	LOC_CALCS
		SimFlag=0;
		rseed=999;
		int on,opt;
		//La linea siguiente evalua si la opcion de linea de comando "-SimFlag" existe
		//Si esta existe usa el numero aleatorio para la simulacion
		if((on=option_match(ad_comm::argc,ad_comm::argv,"-sim",opt))>-1)
		{
			SimFlag=1;
			rseed=atoi(ad_comm::argv[on+1]);
			cout<<"======================\n"<<endl;
			cout<<"      ** Simulacion Poblacional **"<<endl;
			cout<<"======================"<<endl;
		}
	END_CALCS

	//Lee los objetos del archivo de datos desde initprefix
	//Debug
	init_int debug
	//Time indicators
	init_int nmo  //number of months (time index)
	init_int nyr //number of years
	init_int styr //start year
	init_int endyr //end year
	!! cout<<"styr\t"<<styr<<endl;
	!! cout<<"endyr\t"<<endyr<<endl;
	!! cout<<"nmo\t"<<nmo<<endl;

	//Read time-series data from Tindex = 1 to nmo
	init_matrix s_TSdata(1,nmo,1,31);
	//indices and dependent
	vector Tindex(1,nmo);   //Indicador de tiempo secuencial
	vector Month(1,nmo);    //Meses
	vector Years(1,nmo);     //anhos
	vector Effort(1,nmo);
	vector s_obs_catch(1,nmo);   //captura de sardina
	vector s_cpue(1,nmo);
	vector s_obs_Brec(1,nmo); //RECLAS acoustic biomass survey
	vector s_obs_Bpel(1,nmo); //PELACES acoustic biomass survey
	vector s_obs_SSB(1,nmo);  //DEPM spawning biomass
	vector s_aRLP(1,nmo);    //seasonal length-weight parameters
	vector s_bRLP(1,nmo);
	vector s_tempMat(1,nmo);  //seasonal maturity

        //Traspasa datos a vectores 
	LOC_CALCS
		//sardina
		Tindex = column(s_TSdata,1);  //Time index
		Years =  column(s_TSdata,2);   //Years
		Month = column(s_TSdata,3);   //Months
		Effort = column(s_TSdata,8);
		s_obs_catch = column(s_TSdata,6);   //total monthly fishery catch biomass
		s_cpue = column(s_TSdata,10);    //Catch per unit effort
		s_obs_Brec=column(s_TSdata,12);  //Obs acoustic biomass Reclas survey
		s_obs_Bpel=column(s_TSdata,20); //Obs acoustic biomass Pelaces survey
		s_obs_SSB=column(s_TSdata,24);  //DEPM biomass
		s_aRLP=column(s_TSdata,28);     //intercept length-weight relationship
		s_bRLP=column(s_TSdata,29);     //slope LWR
		s_tempMat=column(s_TSdata,26);   //seasonal maturity
	END_CALCS

	//Length indicators
 	init_number s_stlen    //first length class
 	init_number s_steplen  //size step for length classes
 	init_number s_nlen     //number of length classes

	//Read fishery length composition
	init_int s_nobs_fishlen;
	init_ivector s_mo_fishlen(1,s_nobs_fishlen);
	init_matrix s_obs_p_fishlen(1,s_nobs_fishlen,1,s_nlen);
	//init_matrix s_obs_p_fishlen(1,nmo,1,s_nlen)

	//!!if (debug==1){char mark; cout<<"MoFishlen= "<<s_mo_fishlen<<endl; cin>>mark; }  

	//Read Reclas length composition
	init_int s_nobs_reclaslen;
	init_ivector s_mo_reclaslen(1,s_nobs_reclaslen);
	init_matrix s_obs_p_reclaslen(1,s_nobs_reclaslen,1,s_nlen)

	//Read Pelaces length composition
	init_int s_nobs_pelacelen;
	init_ivector s_mo_pelacelen(1,s_nobs_pelacelen);
	init_matrix s_obs_p_pelacelen(1,s_nobs_pelacelen,1,s_nlen)

	vector s_len(1,s_nlen) //marca de clase de longitud

	//BIOLOGY
	init_int s_nages; //number of ages (yearly)
  	init_number s_M;        //Natural mortality rate annual
  	init_number s_Linf;
  	init_number s_Kappa;
  	init_number s_t0;

  	//init_number len1; //First length at age 1
  	init_number age1;				//recruitment age (3 mo for october, 4 mo for november...
  	init_vector s_CVageParms(1,2);	     //Intra-cohort parameters for CV of length at age
  	init_vector s_mat(1,2);	      		     //maturity parameters
  	//init_vector s_fec(1,2);     		     //batch fecundity parameters
	init_int nco; 				     	//number of cohorts per year during the time window
	!! cout<<"NCohorts\t"<<nco<<endl;
	init_ivector mes_rec(1,nco);		//initial month of recruitment cohorts
	init_ivector mes_surv(1,nco);		//last month of surviving cohorts
	init_ivector minageco(1,nco);		//minima edad por cohortes

	!!if (debug==1){char mark; cout<<"minageco= "<<minageco<<endl; cin>>mark; }  

	//some additional
	matrix ages(1,12,1,s_nages);
	number s_len1;
	int j;
	int i;
	int k;
	int l;
	number pi;
	int itmp;
	int surv; //index of survival month of cohort
	int rec;  //index of initial month of cohorts
	imatrix mes(1,12,1,nyr)
	matrix s_ageco(1,nmo,1,nco);  //ages per cohort
	matrix s_lenco(1,nmo,1,nco);   //length at age per cohort
	matrix s_weico(1,nmo,1,nco);	    //weight at age per cohort
	matrix s_matco(1,nmo,1,nco);  //maturity at age per cohort
	matrix ones(1,nmo,1,nco);  //units
	vector tspawn(1,nyr);      //month peak spawning
	vector s_w_age(1,s_nages)
	vector s_mat_age(1,s_nages)

        number s_rho;
	number s_rho_m
	matrix s_sd_age(1,12,1,s_nages)
	matrix s_cv_age(1,12,1,s_nages)
	number diff
	3darray s_Agesize(1,12,1,s_nages,1,s_nlen)
	matrix s_len_age(1,12,1,s_nages)
	number s_Mm; //monthly natural mortality
	number s_nagesm

	LOC_CALCS
		s_Mm=s_M/12;    //Fixed nat mort per month
		s_nagesm=s_nages*12;  //ages months
		//Compute mean length for length-classes
		s_len(1) = s_stlen;
		for(j=1;j<s_nlen;j++)
		{
			s_len(j+1)=s_len(j)+s_steplen;
		}

		//time index for peak spawning biomass (=mes 11=agosto because 1=october)
		//to be used in stock-recruitment
		tspawn(1)=11;
		for(i=1;i<nyr;i++)
		{
			tspawn(i+1)=tspawn(i)+12;
		}
	END_CALCS
	!!if (debug==1){char mark; cout<<"slen"<<s_len<<endl; cin>>mark; }  

	!!if (debug==1){char mark; cout<<"tspawn"<<tspawn<<endl; cin>>mark; }  


	LOC_CALCS
		//age, length-at-age, w-at-age, maturity at age
		for(i=1; i<=nco;i++)
		{
			surv=mes_surv(i);
			rec=mes_rec(i);
			itmp=minageco(i);
			for(j=rec;j<=surv;j++)
			{
				s_ageco(j,i) = itmp;	
				s_lenco(j,i) = s_Linf*(1-exp(-1*s_Kappa*(s_ageco(j,i)-s_t0)/12));
				s_weico(j,i) = s_aRLP(j)*pow(s_lenco(j,i),s_bRLP(j));
				s_matco(j,i) = 1./(1+exp(-(s_mat(1)+s_mat(2)*s_lenco(j,i))));
				ones(j,i)=1.;
				itmp=itmp+1;
			}
		}
	END_CALCS


	LOC_CALCS
		//AGE LENGTH STUFF
		s_len1=s_Linf*(1-exp(-1*s_Kappa*(age1-s_t0)));
		s_len_age(1,1)=s_len1;
		ages(1,1)=age1;
		s_rho=mfexp(-1.*s_Kappa);
		//ages
		//First month
		for(j=1;j<s_nages;j++)
			{
					s_len_age(1,j+1)=s_Linf*(1.-s_rho)+s_rho*s_len_age(1,j);
					ages(1,j+1)=ages(1,j)+1;
			}
		//Now for the rest 11 months
		s_rho_m=mfexp(-1.*s_Kappa/12);
		for(j=1;j<=s_nages;j++)
			{
			for(i=1;i<12;i++)
				 {
					s_len_age(i+1,j)=s_Linf*(1.-s_rho_m)+s_rho_m*s_len_age(i,j);
					ages(i+1,j)=ages(i,j)+0.0833;
				}
			}
			//CV of length at age
			for(i=1;i<=12;i++)
				{
				for(j=1;j<=s_nages;j++)
					{
						s_cv_age(i,j)=(s_CVageParms(1)+s_CVageParms(2)*ages(i,j))/100;
						if(s_cv_age(i,j) < 0){s_cv_age(i,j)=0.002324; }
						//s_cv_age(i,j)=0.05;
						s_sd_age(i,j)=s_cv_age(i,j)*s_len_age(i,j);
					}
				}

			//calcula avg-wt-age and maturity for august peak spw s-r stuff age-structured
			for (j=1;j<=s_nages;j++)
			{
				s_w_age(j) = s_aRLP(11)*pow(s_len_age(11,j),s_bRLP(11));
				s_mat_age(j)=s_tempMat(11)/(1+exp(-(s_mat(1)+s_mat(2)*s_len_age(11,j))));
       			}

			//AGE-LENGTH KEY
			     pi = 3.141592653590;
			     for(k=1;k<=12;k++)
				{
				for(i=1;i<=s_nages;i++)
						{
						for(j=1;j<=s_nlen;j++)
							{
								diff = s_len(j)-s_len_age(k,i);
							     s_Agesize(k,i,j)=(1/sqrt(2*pi*(square(s_sd_age(k,i)))))*mfexp(-0.5*square(diff)/(square(s_sd_age(k,i))));
							}
							s_Agesize(k,i)=s_Agesize(k,i)/sum(s_Agesize(k,i));
						}
				}

				for(i=1;i<=nyr;i++)
				{
					for(j=1;j<=12;j++)
					{
						mes(j,i)=(i-1)*12+j;
					}
				}

					 // cout << mes << endl;
	END_CALCS

	!!if (debug==1){char mark; cout<<"watage= "<<s_w_age<<endl; cin>>mark; }  
	!!if (debug==1){char mark; cout<<"matage= "<<s_mat_age<<endl; cin>>mark; }  

	//End of data file
	init_int eof;	
	LOC_CALCS
	  cout<<"eof = "<<eof<<endl;
	  if(eof==999){
		cout<<"\n ** FIN DE LA SECCION DATA ** \n"<<endl;
	  }else{
		cout<<"\n ¡¡¡¡ ERROR AL LEER LOS DATOS!!!!! \n"<<endl; exit(1);
	  }
	END_CALCS

	number fmsy; //Fishing mortality at Fmsy
	number msy; //Maximum sustainable yield
	number bmsy; //Spawning biomass at MSY
	vector size_tau(1,3); //MLE of variance for size composition
	matrix simcatsize(1,nmo,1,s_nlen); //captura por talla para la simulacion

	//============================================
	// Lee los parametros de control de estimacion
	//============================================
	//Leading parameters
	//theta[1]		log_ro, or log_msy     unexploited recruitment
	//theta[2]		steepness(h), or log_fmsy
	//theta[3]		log_m  : mortality rate
	//theta[4]		log_avgrec : average recruitment
	//theta[5]		prho  : proportion of observation variance recruitment
	//theta[6]		kpsr : slope s-r

	!! ad_comm::change_datafile_name(ControlFile);
	init_int npar;
	init_matrix theta_control(1,npar,1,7)
	
	!!if (debug==1){char mark; cout<<"Control= "<<theta_control<<endl; cin>>mark; }  

	vector theta_ival(1,npar);
	vector theta_lb(1,npar);
	vector theta_ub(1,npar);
	vector theta_phz(1,npar);
	vector theta_prior(1,npar);
	LOCAL_CALCS
		theta_ival = column(theta_control,1);  //ival 
		theta_lb = column(theta_control,2);
		theta_ub = column(theta_control,3);
		theta_phz = column(theta_control,4); //phase
		theta_prior = column(theta_control,5); //prior
	END_CALCS

	init_int ngear; //numero de arte con length frequency
	//=========================
	//       Opciones de selectividad
	//=========================
	// type 1 = logistic (2 parms)
	// type 2 = logistica fija y predeterminada cambiando la phase (sel_phz) a negativa
	init_ivector isel_type(1,ngear);   //tipo de selectividad para la pesqueria, reclas y pelaces
	init_vector ahat(1,ngear);           //age-at-50% vulnerability logistic
	init_vector ghat(1,ngear);	       //std at 50%
	init_ivector sel_phz(1,ngear);	//Phase for estimating selectivity parameters.
	!! cout<<isel_type<<endl;
 
	// fixed logistic (no parameters estimated)
	// ensure sel_phz is set to negative value.
	LOCAL_CALCS
		for(i=1;i<=ngear;i++)
		{
			if(isel_type(i)==2) sel_phz(i) = -1;
		}
		//cout<<"Number of estimated selectivity parameters\n"<<isel_npar<<endl;
	END_CALCS

	//Controls for prior on catchability q (including cpue).
	init_int nits;       //numero de surveys (inits = 4, 1: cpue; 2: reclas; 3: pelaces; 4: mph)
	init_ivector q_prior(1,nits);  //tipo de prior usado para q(0: uniform; 1: normal;..)
	init_vector q_mu(1,nits);   //valor de q assumido escala log
	init_vector q_sd(1,nits);   // std para q
	!! cout<<"q Prior\n"<<q_mu<<endl<<q_sd<<endl;

	init_vector cv_surveys(1,nits); //CVs surveys
	init_vector nsizemult(1,ngear) //effective size for multinomial
	number CV_cpue;
	number CV_reclas;
	number CV_pelaces;
	number CV_mph;
	number n_fish;
	number n_reclas;
	number n_pelaces
	LOCAL_CALCS
		CV_cpue = cv_surveys(1);
		CV_reclas = cv_surveys(2);
		CV_pelaces = cv_surveys(3);
		CV_mph = cv_surveys(4);
		n_fish = nsizemult(1);
		n_reclas = nsizemult(2);
		n_pelaces = nsizemult(3);
	END_CALCS

	//Miscellaneous controls
	// 1 -> verbose
	// 2 -> recruitment model (1=beverton-holt, 2=rickers)
	// 3 -> autocorrelation in recruitment
	// 4 -> std in catch first phase
	// 5 -> std in catch in last phase
	// 6 -> mean fishing mortality rate to regularize the solution
	// 7 -> standard deviation of mean F penalty in first phases
	// 8 -> standard deviation of mean F penalty in last phase.
	
	init_vector cntrl(1,8);
	int verbose;

	init_int eofc;
	LOC_CALCS
		verbose = cntrl(1);
		cout<<"cntrl\n"<<cntrl<<endl;
		cout<<"eofc\t"<<eofc<<endl;
		if(eofc==999){
			  cout<<"\n ** FIN CONTROL ** \n"<<endl;
	  	}else{
			  cout<<"\n ¡¡¡¡ ERROR CONTROL!!!!! \n"<<endl; exit(1);
	 	}
	END_CALCS

	int nf;
	
	ivector ilvec(1,6);
	!!ilvec=3;			//number of gear for selectivity ngear
	!!ilvec(2)=nits;		//number of surveys antes nit
	!!ilvec(3)= 3;         //na_gears number of size-comps
	!!ilvec(4)=1;

PARAMETER_SECTION
	init_number a;
	number b;
	objective_function_value f;
PRELIMINARY_CALCS_SECTION
	//Corre el modelo con parametros conocidos para simular datos
	//nf=0;
	//if(SimFlag)
	// {
	//	initParameters();
	//	calcSelectivities();
	//	calcTotalMortality();
	//	//simulation_model(rseed);
	//}

RUNTIME_SECTION
    maximum_function_evaluations 100,200,500,25000,25000
    convergence_criteria 0.01,0.01,1.e-5,1.e-5

PROCEDURE_SECTION
	b=1.5*a;
	
	//initParameters();

	//calcSelectivities();
	
	//calcTotalMortality();
	
	//calcNumbersAtAge();
	
	//calcFisheryObservations();
		
	//calc_survey_observations();
	
	//calc_stock_recruitment();
	
	//calc_objective_function();

	//sd_depletion=sbt(nyr)/bo;
	
	//if(mceval_phase()) mcmc_output();
	
	f=square(b);


GLOBALS_SECTION
	#undef REPORT
	#define REPORT(object) report << #object "\n" << object << endl;
    
	#undef COUT
	#define COUT(object) cout << #object "\n" << object <<endl;

	#if defined(WIN32) && !defined(__linux__)
		const char* PLATFORM = "Windows";
	#else
		const char* PLATFORM = "Linux";
	#endif

	#include <admodel.h>
	#include <time.h>
	#include <string.h>

	 
	adstring BaseFileName;
	adstring ReportFileName;
	
	adstring stripExtension(adstring fileName)
	{
		/*
		This function strips the file extension
		from the fileName argument and returns
		the file name without the extension.
		*/
		const int length = fileName.size();
		for (int i=length; i>=0; --i)
		{
			if (fileName(i)=='.')
			{
				return fileName(1,i-1);
			}
		}
		return fileName;
	}
